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We describe two new - stochastic-geometrical - methods to obtain reliable velocity 
field statistics from N-body simulations and from any general density and veloc- 
ity fluctuation field sampled at a discrete set of locations. These methods, the 
Voronoi tessellation method and Delaunay tessellation method, are based on the 
use of the Voronoi and Delaunay tessellations of the point distribution defined by 
the locations at which the velocity field is sampled. Adjusting themselves auto- 
matically to the density of sampling points, they represent the optimal estimator 
for volume-averaged quantities. They are therefore particularly suited for checking 
the validity of the predictions of quasi-linear analytical density and velocity field 
perturbation theory through the results of N-body simulations of structure forma- 
tion. We illustrate the subsequent succesfull application of the two methods to 
estimate the bias-independent value of f2 in the N-body simulations on the basis of 
the predictions of perturbation theory for the O-dependence of the moments and 
PDF of the velocity divergence in gravitational instability structure formation sce- 
narios with Gaussian initial conditions. We will also shortly discuss practical and 
complicating issues involved in the obvious extension of the Voronoi and Delaunay 
method to the analysis of observational samples of galaxy peculiar velocities. 



1 Introduction 



The study of the large-scale cosmic velocity field is a very promising and cru- 
cial area for the understanding of structure formation. The cosmic velocity 
field is in particularly interesting because of its close relation to the underlying 
field of mass fluctuations. Indeed, on these large and (quasi)-linear scales the 
acceleration, and therefore the velocity, of any object is expected to have an 
exclusively gravitational origin so that it should be independent of its nature, 
whether it concerns a dark matter particle or a bright galaxy. Moreover, in the 
linear regime the generic gravitational instability scenario of structure forma- 
tion predicts that at every location in the Universe the local velocity is related 
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to the local acceleration, and hence the local mass density fluctuation field, 
through the same universal function of the cosmic density parameter fl (Pee- 
bles 1980), f(fl) oc £! - 6 . Because linear theory provides a good description 
on scales exceeding a few Megaparsec, the use of this straightforward relation 
implies the possibility of a simple inversion of the measured velocity field into 
a field that is directly proportional to the field of local mass density fluctua- 
tions. Such a procedure can then be invoked to infer the value of O, through 
a comparison of the resulting field with the field of mass density fluctuations 
in the same region. 

However, such a determination of Q may be contrived as the estimate 
of the mass density fluctuation field on the basis of the observed galaxy dis- 
tribution may offer a biased view of the underlying mass distribution. By 
lack of a complete and self-consistent physical theory of galaxy formation, the 
commonly adopted approach is to make the simplifying assumption that the 
galaxy density S g and the mass density S are related via a linear bias factor 
b, S g = bS. The comparison between the observed galaxy density fluctuation 
field and the local cosmic velocity field will therefore yield an estimate of the 
ratio {3 — f(Cl)/b w Q°- 6 /b.. However, while numerous studies have yielded 
estimates of /? in the range (3 w 0.5 — 1.2 (see Dekel 1994, Strauss & Willick 
1995, for compilations of results), it has proven very cumbersome to subse- 
quently disentangle the contribution of O and b to the quantity (3. In fact, 
it turns out to be impossible within procedures based on the linearity of the 
analysed velocity field. 

2 Perturbation theory and the Quasi-linear evolution of the veloc- 
ity field 

As long as density and velocity fluctuations are small they grow linear, at a 
global rate irrespective of the value and location of the perturbation. The 
dependence on the mass of the fluctuation only expresses itself through the 
value of the universal cosmic matter density, i.e. through the value of O. This 
specific circumstance will cease to hold as soon as the fluctuations start to 
acquire values in the order of unity and higher. The larger gravitational accel- 
eration exerted by the more massive structures in the density field then starts 
to induce the infall of proportionally ever larger amounts of matter and hence 
to an increase of the growth rate. This leads to a situation in which the growth 
of the density fluctuation is no longer merely dependent on the global value 
of CI, but also on the local value of the density excess or deficit. This in turn 
implies that the situation during the linear regime of the random density and 
velocity perturbations retaining the same statistical properties will no longer 
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Figure 1: The PDF of the velocity divergence 6 for 3 different cosmologies, as a function of 
Q and erg. 

persist. Instead, the nonlinear evolution of the fluctuations leads to an in- 
creasingly asymmetric evolution of the density and velocity fluctuations, with 
overdensities collapsing into compact massive objects whose density excess can 
achieve values exceeding unity by many orders of magnitude, while underden- 
sities expand into empty troughs whose density deficit is naturally limited to 
be no lower than -1.0. 

If the primordial random fluctuation field is characterized by Gaussian 
statistics, one can apply perturbation theory to analytically compute the devel- 
opment of the stochastic properties of the field in the first quasi-linear stages as 
mild nonlinear perturbations develop in the random density and velocity field 
through the action of gravity (see Juszkiewicz & Bouchct 1997 for a review 
and references). For the purpose of determining the value of ft, Bernardeau 
(1994) and Bernardeau et al. (1995) found one of the most straightforward 
and useful results in the context of perturbation theory. In an extensive body 
of work he demonstrated the existence of a simple relation between the higher 
order moments (9 P ) of the divergence 9 of the locally smoothed velocity field 



to its second order moment (9 2 ). The coefficient T p depends on Q, on the shape 
of the power spectrum, the geometry of the window function that has been used 
to filter the velocity and even on the value of the cosmological constant A. For 
the case of a top-hat filtered velocity field, analytical expressions for T p were 
derived that showed them to be strongly dependent on the value of fi, only 
very weakly dependent on A and, very significant within the present context, 
completely independent of a linear bias factor b. As in practical circumstances 
it may still be feasible to get reliable estimates of the lower-order moments, 
this implies that e.g. the relation between the skewness and the second-order 
moment of 9 may be succesfully exploited for the determination of ft through 
the coefficient T 3 , 

T 3 = {9 3 )/{9 2 ) 2 oc IT ' 6 . (2) 

In fact, perturbation theory allows to infer the strongly O-dependent, weakly 
A-dependent, and b- independent expressions for all coefficients T p , and from 
this complete series of coefficients one can construct the complete Probabil- 
ity Distribution Function p{9) (see Bernardeau 1994). An illustration of the 
changing global shape and behaviour of the complete PDF is evidently more 
direct into conveying an impression of the way in which the statistical proper- 
ties of the velocity field depend on the value of Q, and evolve through gravity. 
Such an illustration is provided by the linear-log plot of p{9) for 3 different 
cosmologies. Notice that fl not only influences the overall shape of the PDF 
but also the location of its peak - indicated by 9 maK - as well as that of the 
cutoff at the high positive values of 9. The value of the maximum of 9, at # max 
is directly related to the expansion of the deepest voids in such a Universe. 



3 Volume-averaged quantities 

The perturbation formalism discussed above evidently concerns continuous 
density fields. Practical applications to either N-body simulations or real ob- 
servational samples, on the other hand, yield velocity fields that are sampled 
at a finite number of discrete, non-uniformly distributed, locations. This dis- 
creteness forms a major technical obstacle for the succesfull application of the 
theoretical results. The usual approach is to smooth the discrete velocity field 
by some filter function. Almost without exception the conventional filtering 
schemes concern a mass-weighted velocity field filtered with a filter function 
W M (x,x ) 

J dx v(x) p(x)W M (x, x ) 

v mass( x o) = J • (3) 

/ dxp{x)W M {x,x a ) 
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An example of this is probably one of the most frequently applied class of 
filtering schemes, involving the interpolation of the velocity field values at the 
random sampling locations to those at regular grid locations, weighing the 
contribution by each sampling point by the filter function value. However, the 
presence of the the extra mass- weighting density field factor p would introduce 
considerable technical repercussions, and therefore analytical results within 
perturbation theory have been almost exclusively limited to volume-weighted 
filtered velocity fields v, 



where Wy(x, x ) is the applied weight function. For a succesfull practical 
implementation of the analytical results of perturbation theory it is then of 
crucial importance to have reliable numerical estimators of volume-averaged 
quantities. One specific aspect of such estimators is that from the discrete set 
of sampled velocities they should provide a prescription for the value of the 
velocity field throughout the whole sampling volume, unlike the conventional 
mass-weighted schemes that may restrict themselves to estimates at a finite 
number of positions. 

4 Voronoi and Delaunay Tessellations 

Ideally, the discrete probing and subsequent filtering of the velocity field should 
assign to each location in space the value of the velocity at exactly that location. 
In the infeasible situation of being able to do this all over space this would 
yield the continuous velocity field that one attempts to reconstruct as good as 
possible from the finite number of discretely sampled velocities. In other words, 
one could consider the underlying continuous velocity field as the result of a set 
of an infinite number of sampled points volume filtered with a filter function 
Wy whose filter radius is infinitely small. In the case of a discrete sampling of 
the field this suggests a velocity field reconstruction procedure that defines the 
velocity at every point in space by filtering the corresponding density field 
which can be considered as the sum of delta functions peaking at each sample 
location - with a volume weighted filter that also has a filter radius as small as 
possible. It is then easy to see (see Bernardeau & Van de Weygaert 1996) that 
this defines a velocity field in which the value of the velocity at every location 
in space acquires the value of the velocity at the closest point of the discrete 
velocity sample. 
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Figure 2: Voronoi and Dclaunay tessellations of a 2D set of particles (filled circles). The 
solid lines form the Voronoi tessellation, the dashed lines the Delaunay tessellation. We 
also indicated a normal vector n of the wall separating the points A and B. The grey 
circle represents the area in which one determines the top-hat filtered volume average of the 
velocity gradients. 



The procedure described above implies nothing else than the concept of the 
Voronoi tessellation. Such a tessellation consists of a space-filling network of 
mutually disjunt convex polyhedral cells, the Voronoi polyhedra, each of which 
delimits the part of space that is closer to the defining point in the discrete 
point sample set than to any of the other sample points (see Van de Weygaert 
1991, 1994, for extensive descriptions and references). 

For the application of velocity field analysis, the Voronoi method intro- 
duced by Bernardeau & Van de Weygaert (1996) is based on the assumption 
that the velocity field is uniform within each Voronoi cell of the tessellation, 
such that the velocity throughout each of the Voronoi polyhedra is equal to that 
of the sample point defining that cell. This assumption immediately implies 
that the only non-zero values of the velocity gradients dvi/dxj are localized to 
the (polygonal) Voronoi walls. For the specific case that the window function 
for the volume filtering is a top-hat filter, the subsequent computation of the 
volume averages of the velocity gradients consists of a relatively simple sum 
of the values of those velocity gradients in each of the walls k intersected or 
inside the filter sphere weighted by the surface area Aj~ of the part of the wall 
located within the sphere. Of special interest to us was for instance the lo- 
cally volume-filtered velocity divergence 6>vor , which for a radius R of the filter 
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Figure 3: The PDF of the velocity divergence for various values of Q. The solid lines are the 
theoretical predictions and the dashed lines the predictions for Q = 1 and the same variance. 
The numerical estimates have been obtained using the Delaunay method. 

sphere is computed from 



To illustrate the above we refer to the cartoon representation of the two- 
dimensional Voronoi method in Figure 2. 

While the Voronoi method in general leads to good results and in fact 
was succesfully applied to the results of N-body simulations to confirm the 
predictions of perturbation theory (see Bernardeau & Van de Weygaert 1996), 
it evidently represents an artificial situation of a discontinuous velocity field. 
Moreover, its artificial velocity field implies a few limitations to its application. 
The most significant one of these is that it cannot be applied to filter radii that 
are smaller than the average Voronoi wall distance. Below those scales the 
probability that a randomly placed filter sphere does not contain or intersect 
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any Voronoi wall gets prohibitively high, and therefore it would yield unrealistic 
zero values for the velocity gradients. 

It is therefore best to consider the Voronoi method as a zeroth-order inter- 
polation scheme, and extend the machinery to include a first-order interpola- 
tion scheme. It is the Delaunay method that can be regarded as this elaboration 
and extension towards a multidimensional equivalent of a linear interpolation 
scheme. It is based on the Delaunay tessellation defined by the points in the 
sample. This uniquely defined and space-covering network consists of mutu- 
ally disjunt Delaunay tetrahedra. Each of the Delaunay tetrahedra is defined 
by four nuclei from the point sample that have a circumscribing sphere that 
does not contain any of the other nuclei in its interior. In fact, a host of 
practical, computational and image-processing, applications already employ 
the interpolation and triangulation qualities of Delaunay tessellations, based 
on the realization that they represent the near-optimal multidimensional tri- 
angulation of a discretely sampled space. For better appreciation, we have 
also illustrated the Delaunay method through its two-dimensional equivalent 
in Figure 2, the dashed lines representing the Delaunay tessellation. As they 
are each others dual, the Delaunay and Voronoi tessellation are closely related. 
This close relationship is born out by the fact that the centre of the circum- 
scribing sphere of a Delaunay tetrahedron is a vertex in the corresponding 
Voronoi tessellation. 

The first-order velocity field interpolation scheme defined by our Delaunay 
method (see Bernardeau & Van de Weygaert 1996) is based on the construc- 
tion of the velocity field v(M) at every location M in space through linear 
interpolation between the velocities of the four particles A, B, C and D that 
define the Delaunay cell in which M is situated. v(M) = cya-v(A) +«bv(B) + 
ac'v(C) + a£)v(Z)), where the weights ctj are the barycentric weights of the 
points. The resultant velocity field is one of a field of uniform velocity gradi- 
ents within each Delaunay cell. The value of each of the nine velocity gradient 
components dvi/dvj within each Delaunay tetrahedron, and through them the 
value of the velocity divergence 9, the shear ov,, and even the vorticity u>i, 
follow immediately through solving the 9 linear equations that one obtains 
through evaluation of the three independent velocity differences obtained by 
evaluating the value of the velocity at the four different vertex locations of the 
tetrahedron. Having defined the values of the velocity gradient, and hence of 
the velocity itself, over the entire sample volume, the last step in the Delaunay 
method then consists of determining the volume averaging quantities through 
top-hat filtering with a filter Wth that has a characteristic radius R. Essen- 
tially this consists of determining the weighted average of the value of 9, o~ij 
and/or u>i in a sphere of radius R centered on a location xq, the weights being 
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the volume of the part of the Delaunay tetrahedron that intersects with the 
filter sphere. 

As the velocity gradients in the Delaunay method will nowhere acquire 
artificial zero values, it is a much more robust method. Moreover, it is far less 
memory consuming to store the information on the Delaunay tetrahedra than 
it is to store the complete geometrical information on Voronoi polyhedra, so 
that it can be applied to truly big datasets. The one serious disadvantage in 
its present state is that it is a rather time-consuming method, due to the fact 
that the calculation of the intersection between randomly shaped and located 
tetrahedra and spheres turns out to be anything but a trivial problem. 

By applying the procedures sketched above to a large number of randomly 
located filter spheres of a particular filter radius i?, one obtains a numerically 
determined probability distribution of the velocity gradients in for example 
the outcome of N-body simulations. In this way one can check the analytical 
predictions of perturbation theory, evidently restricted to the early nonlinear 
stages of such simulations. Moreover, as agreement between the analytical 
predictions and the numerical inferences by the Voronoi and Delaunay method 
provides confidence in both, one can apply the tessellation methods to for 
extending the corresponding statistical study to the highly nonlinear structure 
formation stages whose velocity field characteristics cannot be adressed by 
analytical means. 

5 Results, Discussion and Prospects 

The main incentive for developing the Delaunay and Voronoi method is pro- 
vided by the wish to be able to infer a bias-independent value of £1 through 
comparison of the velocity statistics obtained from the discrete point sam- 
ple with those of analytical distributions. In figure 3 we show the PDFs of 
the velocity divergence 9 that were numerically determined by the Delaunay 
method for a range of N-body simulations, each with a different value of fi. 
The solid curve shows the corresponding analytical distribution function p{9), 
for the cosmic epoch with the same dispersion <jg. For contrast, each of the 
four frames also contains the dashed curve for the PDF in an Einstein-de Sitter 
universe with the same value of ae ■ Evidently, the Delaunay method is highly 
succesfull in reproducing the correct statistical distribution, and via the rela- 
tions between the moments of the PDF we indeed obtain very good estimates 
of Q. 

While figure 3 illustrates the potential power of the tessellation meth- 
ods, we are obviously motivated to apply them to more practical situations 
and hence more cumbersome cases where selection and sampling effects and 
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sampling errors are of crucial influence. In particular we hope to be able to 
develop a formalism capable of dealing with observational catalogues of pecu- 
liar velocities of galaxies. In previous work (Bernardeau & Van de Weygaert 
1996, Bernardeau et al. 1997) we already adressed the issue of diluted samples. 
In those cases both methods yielded encouraging results. However, the true 
world will present problems ranging from the fact that one can measure galaxy 
velocities only along the line of sight to complicated selection effects like differ- 
ential Malmquist bias (see e.g. Bcrtschinger et al. 1990, Dekel, Bertschinger 
& Faber 1990). Work on these issues is in progress, but they obviously provide 
a considerable complication. 

Acknowledgments 

We wish to thank Eric Hivon and Francois Bouchet, our collaborators in part 
of the project, and Volker Miillcr for the excellent organization of this Potsdam 
meeting. 

References 

1. F. Bernardeau, ApJ 292, 1 (1992) 

2. F. Bernardeau, R. Juszkiewicz, A. Dekel and F. Bouchet, MNRAS 274, 
20 (1995) 

3. F. Bernardeau and R. van de Weygaert, MNRAS 279, 693 (1996) 

4. F. Bernardeau, R. van de Weygaert, E. Hivon and F. Bouchet, MNRAS, 
1997, in press 

5. E. Bertschinger, A. Dekel, S.M. Faber, A. Dressier and D. Burstein, ApJ 
364, 370 (1990) 

6. A. Dekel, ARAA 32, 371 (1994) 

7. A. Dekel, E. Bertschinger and S.M. Faber, ApJ 364, 349 (1990) 

8. R. Juszkiewicz and F. Bouchet, manuscript, 1997 

9. M.A. Strauss and J.A. Willick, Physics Rep. 261, 271 (1995) 

10. R. van de Weygaert, Ph.D. thesis, Leiden University, 1991 

11. R. van de Weygaert, A&A 283, 361 (1994) 



10 



